# estimate a MA(1)
# y = mu + epsilon_t+theta*epsilon_t-1

# set model: mu=1.2,sigma=1.5,theta = 0.2
rm(list = ls())
e1 <- rnorm(101,sd=1.5)
library(xts)
library(magrittr)
ModelData <- xts(e1,order.by = as.Date(1:101))
names(ModelData) <- 'e1'
ModelData$e2 <- lag(ModelData$e1)
ModelData$y <- 1.2 + ModelData$e1+ 0.2*ModelData$e2

y <- ModelData$y[-1] %>% as.numeric()
# estimation
fn <- function(x){
  # x is parameter vector, include x[1]=mu, x[2]=sigma and x[3]=theta
  epsilon <- numeric(length(y))
  epsilon[1] <- 0
  for (i in 2:length(y)) {
    epsilon[i] <- y[i]-x[1]-x[3]*epsilon[i-1]
  }
  
  L <- length(y)*log(x[2])+ sum(epsilon^2)/(2*x[2]^2) # 似然函数极大化，donlp2默认最小化，加负号
}

library(Rdonlp2)
mle <- donlp2(c(0.1,0.9,0.1),fn)
mle$par
